%% Modify FRAP fitting for each curve
close all;

k = 11; % Select which to correct
Nmean = 10; % Number of final frames to use for final fluorescence value

% Set new FRAP_struct if k=1
if k==1
    FRAP_struct_corrected = FRAP_struct;
end



%% Obtain relevant data vectors from FRAP_struct
time_recovery = FRAP_struct(k).time_recovery; % Recovery time vector (in seconds)
Fcap = FRAP_struct(k).cap_fluor_recovery_corrected; % Fluorescence of cap region (AU)
Fcap_eq = mean(Fcap(end-Nmean:end)); % Equilibrated/final fluorescence


Fcap_norm = (Fcap - Fcap(1))/(Fcap_eq - Fcap(1));

%% Perform fitting
Fcap_fitting = 1 - Fcap_norm;

[f1,f2] = fit(time_recovery,Fcap_fitting,'exp1');

% Compute error in fit coefficients
alpha = 0.95; % 95% confidence interval
ci = confint(f1,alpha); % Obtain intervals of confidence (95%)
df = length(time_recovery); % Degrees of freedom in T-test
t = tinv((1+alpha)/2, df); % Use student's T test distribution
se = (ci(2,:)-ci(1,:)) ./ (2*t); % Standard error of T-test

% Now perform estimation using theoretical Gaussian spot for bleaching
% (D = r^2/4*tauD); r = spot radius, tauD = characteristic recovery
% time
tauD = abs(1./f1.b); % Characteristic recovery time (min)
tauD_err = abs(se(2)/f1.b^2);
Deff = FRAP_struct(k).cap_diameter*FRAP_struct(k).cap_length/(4*tauD); % Effective diffusion coefficient (um^2/s)
Deff_err = tauD_err/tauD^2;

% Convert Deff to viscosity of membrane using brownian diffusion
kB = 1.38E-23; % Boltzmann constant (J/K)
T = 310; % Temperature (K)
R = 2E-9; % Radius of probe molecule in meters (NOT SURE HOW GOOD THIS ESTIMATE IS)

% Calculate effective viscosity based on Brownian diffusion model,
% where D = kB*T/(6*pi*R*eta)
eta_factor = kB*T/(6*pi*R);
SI_conversion_Deff = (1E-12); % To convert Deff from um^2/s to m^2/s
eta = eta_factor/(Deff*SI_conversion_Deff)*10; % Viscosity in Poise (P)
eta_err = Deff_err/Deff^2;
    
%% Plot data and fitting
figure;
plot(time_recovery, Fcap_norm, 'k');
hold on;
plot(time_recovery, 1 - f1(time_recovery), 'r');
xlabel('Time (s)');
ylabel('FRAP recovery');

%% Save new parameters to FRAP_struct_new
FRAP_struct_corrected(k).Deff = Deff;
FRAP_struct_corrected(k).Deff_err = Deff_err;
FRAP_struct_corrected(k).eta = eta;
FRAP_struct_corrected(k).eta_err = eta_err;
FRAP_struct_corrected(k).cap_final_fluor = Fcap_eq;


